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(57) ABSTRACT 

An apparatus and method of analysis for three-dimensional 
(3D) physical phenomena. The physical phenomena may 
include any varying 3D phenomena such as time varying 
polar ice flows. A repesentation of the 3D phenomena is 
passed through a Hilbert transform to convert the data into 
complex form. A spatial variable is separated from the 
complex representation by producing a time based covari- 
ance matrix. The temporal parts of the principal components 
are produced by applying Singular Value Decomposition 
(SVD). Based on the rapidity with which the eigenvalues 
decay, the first 3-10 complex principal components (CPC) 
are selected for Empirical Mode Decomposition into intrin- 
sic modes. The intrinsic modes produced are filtered in order 
to reconstruct the spatial part of the CPC. Finally, a filtered 
time series may be reconstructed from the first 3-10 filtered 
complex principal components. 

46 Claims, 21 Drawing Sheets 
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THREE DIMENSIONAL EMPIRICAL MODE 
DECOMPOSITION ANALYSIS APPARATUS, 
METHOD AND ARTICLE MANUFACTURE 

CROSS REFERENCES TO RELATED 
APPLICATIONS 

This application is a con of application Ser. No. 09/150, 
671 filed Sep. 10, 1998, now U.S. Pat. No. 6,311,130 which 
is a continuation-in-part of application Ser. No. 08/872,586 
filed Jun. 10, 1997, now issued as U.S. Pat. No. 5,983,162, 
which itself claims priority under 35 U.S.C. §11 9(e) to U.S. 
Provisional application Ser. No. 60/023,411 filed on Aug. 
14, 1996 and Ser. No. 60/023,822 filed on Aug. 12, 1996. 

ORIGIN OF INVENTION 

The inventor of the invention described herein is an 
employee of the United States Government. Therefore, the 
invention may be manufactured and used by or for the 
Government for governmental purposes without the pay- 
ment of any royalties thereon or therefor. 

COPYRIGHT NOTIFICATION 

Portions of this patent application contain materials that 
are subject to copyright protection. The copyright owner has 
no objection to the facsimile reproduction by anyone of the 
patent document or the patent disclosure, as it appears in the 
Patent and Trademark Office patent file or records, but 
otherwise reserves all copyright rights whatsoever. 

BACKGROUND OF THE INVENTION 

1. Technical Field of the Invention 

This invention generally relates to a data analysis method, 
apparatus and article of manufacture and more particularly 
to apparatus, article of manufacture and analysis method for 
analyzing three dimensional data varying with respect to a 
forth independent dimension such as three dimensional time 
varying data. 

Although the present invention finds utility in processing 
time varying three dimensional data, it is to be understood 
that any varying n dimensional data (where n^3) represen- 
tative of a real world phenomenon such as data representa- 
tive of a physical process including electrical, mechanical, 
biological, chemical, optical, geophysical or other process 
(es) may be analyzed and thereby more fully understood by 
applying the invention thereto. The real world 
n-dimensional data to which the invention finds utility 
include a wide variety of real world phenomena such as the 
behavior of a stock market, population growth, traffic flow, 
etc. Furthermore, the term “real world n-dimensional data” 
also includes “physical data” representative of physical 
processes such as the electrical, mechanical, biological, 
chemical, optical, geophysical process(es) mentioned above. 

Although the invention is not limited to a particular type 
of signal processing and includes the full range of real world 
data representative of processes or phenomena or combina- 
tions thereof, it is most useful when such real world 
n-dimensional data are nonlinear and non-stationary. 

2. Description of Related Art 

In the parent application, several examples of data from 
geophysical data signals representative of earthquakes, 
ocean waves, tsunamis, ocean surface elevation and wind 
were processed to show the invention’s wide utility to a 
broad variety of signal and data types. The techniques 
disclosed therein and elaborated upon herein represent major 
advances in physical data processing. 
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Previously, analyzing data, particularly those having non- 
linear and/or nonstationary properties, was a difficult prob- 
lem confronting many industries. These industries have 
harnessed various computer implemented methods to pro- 
5 cess data measured or otherwise taken from various pro- 
cesses such as electrical, mechanical, optical, biological, and 
chemical processes. Unfortunately, previous methods have 
not yielded results which are physically meaningful. 

Among the difficulties found in conventional systems is 
10 that representing physical processes with physical signals 
may present one or more of the following problems: 

(a) The total data span is too short; 

(b) The data are nonstationary; and 

(c) The data represent nonlinear processes. 

15 Although problems (a)-(c) are separate issues, the first 
two problems are related because a data section shorter than 
the longest time scale of a stationary process can appear to 
be nonstationary. Because many physical events are 
transient, the data representative of those events are nonsta- 
20 tionary. For example, a transient event such as an earthquake 
will produce nonstationary data when measured. 
Nevertheless, the nonstationary character of such data is 
ignored or the effects assumed to be negligible. This 
assumption may lead to inaccurate results and incorrect 
25 interpretation of the underlying physics as explained below. 

A variety of techniques have been applied to nonlinear, 
nonstationary physical signals. For example, many com- 
puter implemented methods apply Fourier spectral analysis 
to examine the energy-frequency distribution of such sig- 
30 nals. 

Although the Fourier transform that is applied by these 
computer implemented methods is valid under extremely 
general conditions, there are some crucial restrictions: the 
system must be linear, and the data must be strictly periodic 
35 or stationary. If these conditions are not met, then the 
resulting spectrum will not make sense physically. 

A common technique for meeting the linearity condition 
is to approximate the physical phenomena with at least one 
linear system. Although linear approximation is an adequate 
40 solution for some applications, many physical phenomena 
are highly nonlinear and do not admit a reasonably accurate 
linear approximation. 

Furthermore, imperfect probes/sensors and numerical 
schemes may contaminate data representative of the phe- 
45 nomenon. For example, the interactions of imperfect probes 
with a perfect linear system can make the final data nonlin- 
ear. 

Many recorded physical signals are of finite duration, 
nonstationary, and nonlinear because they are derived from 
50 physical processes that are nonlinear either intrinsically or 
through interactions with imperfect probes or numerical 
schemes. Under these conditions, computer implemented 
methods which apply Fourier spectral analysis are of limited 
use. For lack of alternatives, however, such methods still 
55 apply Fourier spectral analysis to process such data. 

In summary, the indiscriminate use of Fourier spectral 
analysis in these methods and the adoption of the stationarity 
and linearity assumptions may give inaccurate results some 
of which are described below. 

60 First, the Fourier spectrum defines uniform harmonic 
components globally. Therefore, the Fourier spectrum needs 
many additional harmonic components to simulate nonsta- 
tionary data that are nonuniform globally. As a result, energy 
is spread over a wide frequency range. 

65 For example, using a delta function to represent the flash 
of light from a lightning bolt will give a phase-locked wide 
white Fourier spectrum. Here, many Fourier components are 
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added to simulate the nonstationary nature of the data in the 
time domain, but their existence diverts energy to a much 
wider frequency domain. Constrained by the conservation of 
energy principle, these spurious harmonics and the wide 
frequency spectrum cannot faithfully represent the true 
energy density of the lighting in the frequency and time 
space. 

More seriously, the Fourier representation also requires 
the existence of negative light intensity so that the compo- 
nents can cancel out one another to give the final delta 
function representing the lightning. Thus, the Fourier com- 
ponents might make mathematical sense, but they often do 
not make physical sense when applied. 

Although no physical process can be represented exactly 
by a delta function, some physical data such as the near field 
strong earthquake energy signals are of extremely short 
duration. Such earthquake energy signals almost approach a 
delta function, and they always give artificially wide Fourier 
spectra. 

Second, Fourier spectral analysis uses a linear superpo- 
sition of trigonometric functions to represent the data. 
Therefore, additional harmonic components are required to 
simulate deformed wave profiles. Such deformations, as will 
be shown later, are the direct consequence of nonlinear 
effects. Whenever the form of the data deviates from a pure 
sine or cosine function, the Fourier spectrum will contain 
harmonics. 

Furthermore, both nonstationarity and nonlinearity can 
induce spurious harmonic components that cause unwanted 
energy spreading and artificial frequency smearing in the 
Fourier spectrum. The consequence is incorrect interpreta- 
tion of physical phenomena due to the misleading energy- 
frequency distribution for nonlinear and nonstationary data 
representing the physical phenomenon. 

According to the above background, the state of the art 
does not provide a useful computer implemented tool for 
analyzing nonlinear, nonstationary physical signals. Geo- 
physical signals provide a good example of a class of signals 
in which this invention is applicable. Great grandparent 
application Ser. No. 08/872,586 filed on Jun. 10, 1997, now 
issued as U.S. Pat. No. 5,983,162 illustrates several types of 
nonlinear, nonstationary geophysical signals which are very 
difficult to analyze with traditional computer implemented 
techniques including earthquake signals, water wave 
signals, tsunami signals, ocean altitude and ocean circulation 
signals. 

This application extends the technique of the parent, 
grandparent and great-grandparent applications to the pro- 
cessing of three-dimensional data and signals representative 
thereof. Three-dimensional data representing, for example, 
time varying physical phenomena are an increasing subject 
of various processing techniques. In fact, the above- 
described prior art techniques such as Fourier analysis are 
routinely applied to process such three-dimensional data. 

Three-dimensional physical phenomena often are nonlin- 
ear and/or nonstationary. Therefore, like the one- 
dimensional and two-dimensional data processing tech- 
niques described above, the conventional processing 
techniques are simply inadequate to process such three- 
dimensional data. 

Moreover, conventional data analysis methods are uti- 
lized to separate the various scales contained in the three- 
dimensional data. For example, in antarctic ice flow and 
antarctic ice protuberance evaluation processing, conven- 
tional spectral analysis was applied toward an objective 
analysis of the information contents. However, where the 
data contained spatial and temporal components, the pres- 
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ence or absence of spatial propagation was not determinable 
and led to erroneous results. 

SUMMARY OF THE INVENTION 

5 It is a purpose of the present invention is to solve the 
above-mentioned problems in conventional signal analysis 
techniques; 

It is another purpose of the present invention to analyze 
three-dimensional phenomena varying with respect to a 
10 fourth independent dimension; 

It is yet another purpose of the present invention to 
perform spectral analysis on non-stationary three- 
dimensional geophysical phenomenon data, such as a time 
sequence of gridded geophysical data maps of a particular 
15 geographical location. 

The present invention permits valid intercomparison of 
various data sequences on an intrinsic mode-by-mode basis 
to enhance greatly the ability to establish teleconnections 
between various geophysical variables. The present inven- 
tion is an apparatus, article of manufacture and method of 
analysis wherein a physical phenomenon is analyzed by: (1) 
removing a linear trend line to keep subsequent matrix 
operations within machine capabilities; (2) passing the data 
25 through an Hilbert transform to convert the data into com- 
plex form; (3) removing a spatial variable by producing a 
covariance matrix in time; (4) producing the temporal parts 
of the principal components by applying Singular Value 
Decomposition (SVD); (5) based on the rapidity with which 
the eigenvalues decay, select the first 3-10 principal com- 
ponents (PCs) for Empirical Mode Decomposition (EMD) 
into intrinsic modes; (6) selecting the intrinsic modes pro- 
duced in order to produce either a high-pass, low-pass, or 
bandpass filter; (7) reconstructing a spatial part of the PCs 
35 (multiplied by the appropriate eigenvalues) by multiplying 
the unfiltered temporal part into the original data; and, as a 
final check, (8) reconstructing a filtered time series from the 
first 3-10 filtered complex PCs. 

Computer implemented Empirical Mode Decomposition 
40 method decomposes principal components representative of 
a physical phenomenon into Intrinsic Mode Functions 
(IMFs) that are indicative of intrinsic oscillatory modes in 
the physical phenomenon. 

Contrary to almost all the previous methods, this new 
45 computer implemented method is intuitive, direct, a 
posteriori, and adaptive, with the basis of the decomposition 
based on and derived from the physical signal. The bases so 
derived have no close analytic expressions, and they can 
only be numerically approximated in a specially pro- 
50 grammed computer by utilizing the inventive methods dis- 
closed herein. 

More specifically, using the general method of the present 
invention the physical signal is analyzed without suffering 
the problems associated with computer implemented Fourier 
55 analysis. Namely, the present invention avoids inaccurate 
interpretation of the underlying physics caused in part by 
energy spreading and frequency smearing in the Fourier 
spectrum. 

Further scope of applicability of the present invention will 
60 become apparent from the detailed description given here- 
inafter. However, it should be understood that the detailed 
description and specific examples, while indicating pre- 
ferred embodiments of the invention, are given by way of 
illustration only, since various changes and modifications 
65 within the spirit and scope of the invention will become 
apparent to those skilled in the art from this detailed descrip- 
tion. Furthermore, all the mathematic expressions are used 
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as a short hand to express the inventive ideas clearly and are 
not limitative of the claimed invention. 

BRIEF DESCRIPTION OF DRAWINGS 

The present invention will become more fully understood 
from the detailed description given hereinbelow and the 
accompanying drawings which are given by way of illus- 
tration only, and thus are not limitative of the present 
invention, and wherein: 

FIG. 1 is a high-level block diagram of a computer system 
which may be programmed with the inventive with the result 
being a special purpose computer; 

FIG. 2 is a high-level flowchart describing the overall 
inventive method which may be implemented on the com- 
puter system shown in FIG. 1; 

FIG. 3 (a) is a high-level flowchart describing the Sifting 
Process which may be implemented on the computer system 
shown in FIG. 1; 

FIG. 3(b) is a continuation of the high-level flowchart in 
FIG. 3(a) describing the Sifting Process which may be 
implemented on the computer system shown in FIG. 1; 

FIG. 4(a) is a graph of data with intermittency for 
illustrating an optional intermittency test of the invention; 

FIG. 4(b) is a graph illustrating the upper envelope, lower 
envelope, envelope mean and original data which are uti- 
lized to explain the computer implemented Empirical Mode 
Decomposition method of the invention; 

FIGS. 4(c) is a graph of the first intrinsic mode function 
when the Sifting Process is applied to the data of FIG. 4(a) 
without applying the intermittency test option; 

FIGS. 4 (d)-(f) are graphs of the first, second, and third 
intrinsic mode functions when the Sifting Process is applied 
to the data of FIG. 4(a) which applies the intermittency test 
option. 

FIGS. 4(g) shows a graph of the first component signal h ± 
generated by the Sifting Process of the invention; 

FIG. 4(h) shows wind speed data in the form of a graph 
plotting wind speed as a function of time for explaining the 
computer implemented Empirical Mode Decomposition 
method of the invention; 

FIG. 4(z) is a graph of the first intrinsic mode function 
component which is generated by the Sifting Process of the 
invention; 

FIG. 5 shows Antarctic sea ice extents in 1° meridional 
sectors on Sep. 1, 1996; 

FIGS. 6(a)-(b) show real spatial and temporal parts of the 
first 3 Complex Principal Components of the time series of 
the Antarctic sea ice extents in 1° meridional sectors; 

FIGS. 7 (a)-(b) show Eigenvalue ratios of the first 20 
complex principal components (CPCs) and the Percent of 
signal in the first N principal components (PCs), obtained 
from the ratio of the sum of the first N eigenvalues to their 
total sum; 

FIG. 8 shows intrinsic modes of the real temporal part of 
the first CPC; 

FIGS. 9 (a)-(b) are a contour plot of the low-pass (lp) 
filtered data of the sectored Antarctic sea ice extents recon- 
structed from the sum of the intrinsic modes with periods 
longer than 3 years of the first 3 Principal Components and 
a three-dimensional plot of the same lp -filtered data; 

FIGS. 10(a)-(Z?) show a contour plot of the intrinsic 
seasonal mode of the sectored Antarctic sea ice extents 
reconstructed from the quasiquadrennial (QQ) parts of the 
first 3 PCs and a three-dimensional plot of the same physical 
phenomenon data. 
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DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

Computer for Implementing the Inventive Method 

Turning now to the drawings, a computer system 50 and 
5 computer 52 suitable for programming with the inventive 
method is diagrammatically shown in the block diagram of 
FIG. 1. To allow human interaction with the computer 52, 
the computer system includes a keyboard 54 and mouse 56. 
The computer, programmed with the inventive method, is 
10 analogous to an electrical filter or a mechanical sieve: it 
separates digital data into series of IMF’s according to their 
time scales in a manner analogous to electrically filtering 
harmonics or a mechanical sieve which separates aggregated 
sand particles according to their physical size. 

Because the invention is applied to analyze physical 
signals, the computer system 100 also includes an input 
device 58, and/or an imaging device 60 which are used to 
collect image information on a physical phenomenon and 
generate physical data representative thereof. Input device 
20 58 may include an antenna connected to a radio receiver for 
receiving radio transmitted data, for example, from a remote 
imaging station such as the National Aeronautics and Space 
Administration (NASA) Nimbus 7 satellite. 

To output the results of the computer implemented 
25 method, the computer system 50 also includes a display 62 
such as a cathode ray tube or flat panel display, a printer 64 
and an output device 66. Each of these outputs (62, 64, 66) 
should have the capability to generate or otherwise handle 
color outputs because, for example, the Hilbert Spectrum 
30 may be in color. 

The generalized output device 66 may also include a 
network connection to connect the computer 50 to a local or 
wide area network. In this way, the physical signal may be 
inputted from the network. Thus, all outputs can be sent to 
35 another location via such a network connection. 

Furthermore, the computer system 50 also includes a 
mass storage device 68. The mass storage device 68 may be 
a hard disk, floppy disc, optical disc, etc. The mass storage 
device 68 may be used to store a computer program which 
40 performs the invention when loaded into the computer 52. 
As an alternative, the input device 58 may be a network 
connection or off-line storage which supplies the computer 
program to the computer 52. 

More particularly, the computer program embodiment of 
45 the invention may be loaded from the mass storage device 68 
into the internal memory 70 of the computer 52. The result 
is that the general purpose computer 52 is transformed into 
a special purpose machine that implements the invention. 

Even more particularly, each step of inventive method 
50 will transform at least a portion of the general purpose 
computer 52 into a special purpose computer module imple- 
menting that step. For example, when a sifting step is 
implemented on the computer 52, the result is a computer 
implemented sifting apparatus (sifter) that performs the 
55 sifting functions of sifting step. 

Other embodiments of the invention include firmware 
embodiments and hardware embodiments wherein the 
inventive method is programmed into firmware (such as 
EPROM, PROM or PLA) or wholly constructed with hard- 
60 ware components. Constructing such firmware and hardware 
embodiments of the invention would be a routine matter to 
one of ordinary skill using known techniques. 

Article of Manufacture 

Still further, the invention disclosed herein may take the 
65 form of an article of manufacture. More specifically, the 
article of manufacture is a computer-usable medium, includ- 
ing a computer-readable program code embodied therein 
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wherein the computer-readable code causes computer 52 to 
execute the inventive method. 

A CD-ROM, a DVD-ROM, high density moveable stor- 
age or a computer diskette such as disc 72 in FIG. 1 are 
examples of such a computer-usable medium. When the disc 
72 is loaded into the mass storage device 72, the computer- 
readable program code stored therein is transferred into the 
computer 52. In this way, the computer 52 may be instructed 
to perform the inventive methods disclosed herein. 

Method for Carrying Out the Invention 

Before describing the computer implemented Three- 
Dimensional (3-D) Empirical Mode Decomposition (3D- 
EMD) analysis method in detail, the definition and physical 
meaning of intrinsic mode functions in general is discussed. 
Intrinsic Mode Function 

An Intrinsic Mode Function (IMF) is a function that 
satisfies the following two conditions: 

(a) in the whole data set, the number of extrema and the 
number of zero -crossings must either be equal or differ 
at most by one, and 

(b) at any point, the mean value of upper envelope defined 
by the maxima and the lower envelope defined by the 
minima is zero. 

The first condition shares some similarity to the tradi- 
tional narrow band requirements for a stationary Gaussian 
process. The second condition is a totally new idea. 
Conceptually, the second condition modifies the classical 
global requirement to a local one. Furthermore, the second 
condition has the desirable result that the instantaneous 
frequency will not have unwanted fluctuations induced by 
asymmetric wave forms. Mathematically, the second condi- 
tion should ideally be “the local mean of the data being 
zero.” For nonstationary data, the “local mean” requires a 
“local time scale” to compute the mean, which is impossible 
to define. Fortunately, the local time scale need not be 
defined to fulfil the second condition, as is discussed below. 

To apply these concepts to physical data, the local mean 
of the signal envelopes are used to force the local symmetry. 
The signal envelopes are defined by the local maxima and 
the local minima. This is an approximation which avoids the 
definition of a local averaging time scale. With the physical 
approach and the approximation adopted here, the inventive 
method does not always guarantee a perfect instantaneous 
frequency under all conditions. Nevertheless, it can be 
shown that, even under the worst conditions, the instanta- 
neous frequency so defined is still consistent with the 
physics of the system being studied and represents the 
system being studied much more accurately than previous 
techniques based on Fourier analysis. 

The term “Intrinsic Mode Function” is adopted because it 
represents the oscillation mode embedded in the data. With 
this definition, the IMF in each cycle, defined by the 
zero-crossings, involves only one mode of oscillation. In 
other words, each IMF represents only one group of oscil- 
lation modes or time scales and no riding waves are allowed 
in an IMF. 

Before presenting the inventive 3D-EMD analysis 
method for decomposing a physical phenomena represented 
as three-dimensional data into IMFs, a qualitative assess- 
ment of the intrinsic oscillatory modes may be roughly 
determined by simply examining the data by eye. From this 
examination, one can immediately identify the different 
scales directly in two ways: the time lapse between the 
successive alternations of local maxima and minima and the 
time lapse between the successive zero -crossings reveals the 
different scales. The interlaced local extrema and zero- 
crossings give us complicated data: one undulation is riding 
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on top of another, and they, in turn, are riding on still other 
undulations, and so on. Each of these undulations defines a 
characteristic scale or oscillation mode that is intrinsic to the 
data: hence, the term “Intrinsic Mode Function” is adopted. 
5 To reduce the data into the needed IMFs, the invention 
utilizes a computer implemented Empirical Mode Decom- 
position Method which is described hereinbelow. 

Overview of Three-Dimensional Signal Processing 

The term “three-dimensional data” is used herein to 
denote the measurement of a physical quantity across three - 
dimensions. For example, three-dimensional data include 
images generated by a multichannel microwave imagers on 
board the NASA Nimbus 7 and the DOD DMSP F 8 , Fll and 
F13 satellites. 

FIG. 2 illustrates the overall inventive method of process- 
ing three-dimensional data including generation of intrinsic 
mode functions (the Sifting Process). To begin the preferred 
processing method, the data are collected for the three - 
20 dimensional physical activity, process or phenomena, sensed 
by an appropriate sensor or imager in step 100. For example 
an image of a physical phenomena may be sensed by a 
multichannel microwave imager 58 to generate a three- 
dimensional representation. 

25 In this example, the three-dimensional data represent 
propagation of sea ice protuberances around Antarctica. The 
data are analyzed using the analysis technique of the present 
invention, which is appropriate for non-stationary time 
series, in particular using empirical decomposition to extract 
30 intrinsic mode functions from the data. Previously, spectral 
analysis was applied to the interior of the Antarctic sea ice 
pack using band-limited discreet Fourier transform of a 
shorter 1978-1987 sea ice time series. This prior approach 
provided phase distributions of two components observed in 
35 the El Nino Southern Oscillation (ENSO) signal, with peri- 
ods of ca. 2.4 and 4.2 years, suggestive of a propagating 
wave in the interior of the pack. 

By contrast, the analysis method of the present invention 
examines sea ice extents, in this example, using Empirical 
40 Mode Decomposition (EMD) in combination with Complex 
Singular Value Decomposition (CSVD). First, CSVD is used 
to separate the data into spatial and temporal components 
and, then, EMD is used to divide the temporal signal into 
distinct modes, representing narrow frequency pass-bands. 
45 It is understood that as used herein complex is intended to 
have its ordinary engineering/scientific meaning, i.e., having 
a real and an imaginary or phase component, e.g., imped- 
ance may include resistance and reactance. Reconstruction 
of each intrinsic mode from its first 3 CPC’s then permits 
50 construction of contour and 3-D plots that reveal the pres- 
ence or absence of spatial propagation of respective nodes. 

So, in step 100 of FIG. 2, the 3D data are collected in a 
sequence of grid maps, A(x,y;t). In step 102 the 3D data are 
reshaped into a 2-D array AA(x*y;t). In the method of 
55 preferred embodiment, the 3D representation is segmented 
into sectors identical at each time sample and extents are 
calculated for each particular sector. The extents are deter- 
mined essentially by summing all of the area within the 
sector above a threshold level (e.g., 15%) for the particular 
60 characteristic being analyzed, e.g., ice concentrations. 
Typically, the area within each sector is segmented into 
smaller identically sized areas of a given granularity and the 
number of smaller areas above the threshold are counted. So, 
for example, representing each of the smaller areas as a pixel 
65 on a display, an area of polar ice may be represented as one 
degree longitudinal sectors of sea ice extents (the sum of 
areas of pixels with ice concentrations at or above 15%). 
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Next, in step 104, the reshaped 3D data may be 
deseasoned by subtracting a model seasonal cycle and a 
trend line, which are both obtained using multiple linear 
regression to model the data. Thus, the deseasoned 
(AA Deseasoned) result may be represented: 

5 

AA Deseasoned = AA - Op + ajt + y] \a 2k cos(27rkt / t ) + +i )Sixi(2jrA:r / r)], 

k= 1 

where a n are coefficients of the model obtained as a result of 
the multiple linear regression and t is the annual period in 
units of time (t). Deseasoning the data may help to contain 
the subsequent matrix operations within the operating capa- 
bilities of the particular machine. 

Continuing in step 106, the deseasoned data are passed 
through an Hilbert transform to preserve spatial phase 
information, A comp/ex =hilb (AA Deseasorled ). Then, in step 
108, the complex Principal Components (PC) are obtained 
using Singular Value Decomposition (SVD), separating the 
spatial component of the Hilbert transform, U, from the 
temporal component, V. Thus, 

A =TT*S*V' 

where, 

U=U(position, PC#); 

S=S(PC#,PC#) (eigenvalues, or weights) and 

V=V(time, PC#). 
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in 1 -degree longitudinal sectors for the period 1978-1996, 
as observed with the multichannel microwave imagers on 
board the NASA Nimbus 7 and the DOD DMSP F8, Fll, 
and FI 3 satellites. After converting the time series into 
5 complex numbers by means of a Hilbert transform, the time 
series of the 360 sectors was decomposed into its complex 
principal components (CPCs), effectively separating the 
spatial and temporal values. Then, the real and imaginary 
parts of the temporal portions of the first 3 CPCs were 
10 decomposed into their intrinsic modes using Empirial Mode 
Decomposition. Each decomposition represents a narrow 
frequency band, resulting in a collection of 3 CPCs for each 
intrinsic mode. Finally, the data were filtered and recon- 
15 structed in two different ways. First, using a low-pass (lp) 
filter on the data, all of the intrinsic modes of each CPC with 
periods longer than 2 years were combined and the result 
designated as lp-filtered. Next, the intrinsic mode of each 
CPC with periods of approximately 4 years was selected and 
20 designated the quasiquadrennial (QQ) modes. The lp-filtered 
time series showed eastward-propagating azimuthal motion 
in the Ross and Weddell Seas, but no clearly circumpolar 
motion. The QQ time series, on the other hand, clearly 
showed eastward-propagating circumpolar waves, but with 
25 occasional retrograde motion to the west. 

Thus the steps described above were implemented in 
MatLab language as follows: 


100 Let the gridded data series by represented by D(x,y,t) 

102 DD(x * y,t)=detrend(D(x*y,t)) /*The spatial coordinates have been 

arranged in a one-dimensional vector*/ 

104 hDD=hilbert(DD) 

106 hDDD=transpose (hDD)*hDD 

108 [ V, SS, V]=svd(hDDD) /^Procedure to produce the temporal parts of 

the PCs, V, and the eigenvalues, S. Because 
hDDD is symmetric, V appears twice on the 
left-hand side */ 

110 VV(time,PC#,mode)=emd(V(time,PC#)) 

/*EMD step 110 applied to the temporal part 
of the first 3-10 PCs.*/ 

112 VVf=transpose(W(time,PC#,selected modes)) 

/* e.g., selected modes are a single mode or, 
all modes with period > one year */ 

114 U(position,selected PCs)*S(selected PCs, selected PCs)=hDD*V(time, selected PCs) 

/* Produces an unfiltered U */ 

116 DD= real(U(position,selected PCs)*S(selected PCs, selected PCs)* 

VV(time, selected PCs, selected PCs) 

/*Reconstructs a filtered time series */ 


In step 110, the temporal parts of the first 10 CPCs, for 50 
example, are decomposed into their respective intrinsic 
modes by means of Empirical Mode Decomposition. The 
prime on V in the above equation for A complex indicates the 
matrix transpose of V. In step 112, the data are reconstructed 
as filtered by selecting a single mode or collection of modes 55 
(the first ten, for example) constituting a band-pass filter (a 
low pass filter in this example): 

A complex , ^ ered =!7(position, 1:10) *5(1:10, 

1 : 10 )’. 

Finally, in step 114, the real part of A complex _ filtered is 
reshaped into a sequence of filtered grid maps and the 
filtered grid map results are displayed in 116. 

The preferred method of the present invention has been 
applied successfully to an Antarctic circumpolar wave in the 65 
perimeter of the ice pack, represented as a time series of the 
sea ice extents (essentially the area within the ice perimeter) 


The Hilbert Spectrum 

Prior to the sifting process of step 110, the Hilbert 
Transform is applied in step 106 to each component to 
preserve spatial phase information A complex . 

Thus, a brief tutorial on the Hilbert transform with empha- 
sis on its physical interpretation can be found in Bendat and 
Piersol, 1986, “Random Data: Analysis and Measurement 
Procedures ” 2nd Ed., John Wiley & Sons, New York, N.Y. 
Essentially, the Hilbert transform is the convolution of X(t) 
with 1/t. This convolution emphasizes the local properties of 
60 X(t). In more formal terms, given a time series X(t), the 
Hilbert Transform Y(t) can be expressed as 


where P indicates the Cauchy principal value. 
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With this definition, X(t) and Y(t) form a complex con- 
jugate pair. This complex conjugate pair Z(t) may be 
expressed as: 

Z(t)=X(t)+iY(t)=a(t)e ie W, 
in which 

a(i)=[X 2 (t)+Y 2 (t)] 1/2 , 

X(t) 

0(l) = arctan— — . 


So, in step 106, the Hilbert transform is applied to the 
deseasoned data (AA Deseasone J) to provide A complex m 
which preserves spatial phase information and contains the 
temporal and spatial information for each of the CPCs. SVD 
is applied to A comp i ex in step 108 to extract the CPCs. EMD 
is applied in step 110 individually to each of these extracted 
CPCs. 

Empirical Mode Decomposition (EMD): The Sifting Pro- 
cess 

FIGS. 3A-B show the steps in the Empirical Mode 
Decomposition of step 110 wherein CPCs are decomposed 
into their respective intrinsic modes. The essence of the 
EMD step 110 is to identify, empirically, the intrinsic 
oscillatory modes by their characteristic time scales in the 
data, and then decompose the data accordingly. The decom- 
position is based on the following general assumptions: 

a. the signal has at least two extrema: one maximum and 
one minimum, and 

b. the characteristic timescale is defined by the time lapse 
between the extrema. 

In other words, the time lapse between successive extrema 
are employed as the definition of the time scale for the 
intrinsic oscillatory mode because it gives a much finer 
resolution of the oscillatory modes and because it can be 
applied to data with non-zero mean (either all positive or all 
negative values, without zero -crossings). A systematic way 
to extract the intrinsic mode functions is the computer 
implemented Empirical Mode Decomposition method or 
Sifting Process which is described as below. 

Optional intermittency tests (201,221) may be introduced 
to alleviate the aliasing effects associated with intermittence 
in the data that can cause mode mixing. These optional 
intermittency tests ( 201 , 221 ) are conducted before con- 
struction of the envelope in steps 210 and 230. Optional 
intermittency test 201 checks the distance between succes- 
sive maxima for whether the distance is within a pre- 
assigned value n times the shortest distance between waves. 
If not, then an intermittency exists and the next step is step 
202. Otherwise, there is no intermittency and the upper 
envelope is constructed in step 210 as further described 
below. 

Similarly optional intermittency test 221 checks the dis- 
tance between successive minima for whether the distance is 
within a pre -assigned value n times the shortest distance 
between waves. If not, then an intermittency exists and the 
next step is step 222. Otherwise, there is no intermittency 
and the upper envelope is constructed in step 230 as further 
described below. 

An example of such intermittency is shown in FIG. 4(a), 
in which small scale waves appear only intermittently. By 
strict application of the Sifting Process, the resulting IMF’s 
are given in FIGS. 4(b), in which two drastically different 
time scales are present in the first IMF component as shown 
in FIG. 4(c). This mixing of modes breaks up the main wave 
train by the small intermittent oscillations. 
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If intermittency tests (201, 222) are employed that utilize 
a preassigned value of n-times the shortest distance between 
waves, the resulting IMFs are shown in FIGS. 4 (d)-(f), in 
which the modes are clearly and physically separated. 
5 Effectively, mode mixing is eliminated by treating the local 
extrema that failed the intermittency test ( 201 , 221 ) as local 
maxima and minima (in steps 202 and 212 ), respectively. 
Therefore, the upper and lower envelope will be identical as 
the original data reference line. 
io As noted above, these intermittency tests 201, 221 and 
subsequent steps 202, 222 are optional. By selecting an 
artificially large value for n in steps 201 and 221 to test for 
intermittency, the test will be effectively bypassed. 
Otherwise, these tests 201, 221 can be bypassed at the initial 
15 selection of the program. 

Next, an upper envelope 20 is constructed for the physical 
signal 10 in step 210. The upper envelope 20 is shown in 
FIG. 4(b) using a dot-dashed line. This upper envelope 20 is 
preferably constructed with a cubic spline that is fitted to the 
20 local maxima. 

In step 220, the local minimum values of the physical 
signal 10 are identified. To complete the envelope, a lower 
envelope 30 is constructed from the local minimum values 
in step 230. The lower envelope 30 is also shown in FIG. 
25 4(b) using a dot-dash line. Like the upper envelope 20, the 
lower envelope 30 is preferably constructed with a cubic 
spline that is fitted to the local minima. 

The upper and lower envelopes 20, 30 should encompass 
all the data within the physical signal 10. From the upper and 
30 lower envelopes 20, 30, an envelope mean 40 is the deter- 
mined in step 240. The envelope mean 40 is the mean value 
of the upper and lower envelopes 20, 30. As can be seen in 
FIG. 4(b), the envelope mean 40 bisects the physical signal 
10 quite well. 

35 Then, the first component signal h-L is generated in step 
250 by subtracting the envelope mean 40 from the physical 
signal 10. This computer implemented step 250 may also be 
expressed as: 

40 Xty-m^h-t, 

where the envelope mean 40 is m 1 and the physical signal is 
X(t). 

FIG. 4(g) shows the first component signal h ± . Ideally, the 
first component signal h a should be an IMF, for the con- 
45 struction of h ± described above seems to have made h 2 
satisfy all the requirements of an IMF. In reality, however, a 
gentle hump that resides on a slope region of the data can 
become an extremum when the reference coordinate is 
changed from the original rectangular coordinate to a cur- 
50 vilinear coordinate. For example, the imperfection of the 
envelopes 20, 30 can be seen by observing the overshoots 
and undershoots at the 4.6 and 4.7 second points in FIG. 
4(b). 

An example of this amplification can be found in the 
55 gentle hump between the 4.5 and 4.6 second range in the 
data in FIG. 4(b). After the first round of sifting, the gentle 
hump becomes a local maximum at the same time location 
in the first component signal h 2 shown in FIG. 4(b). New 
extrema generated in this way actually recover the proper 
60 modes lost in the initial examination. Thus, the Sifting 
Process extracts important information from the signal 
which may be overlooked by traditional techniques. In fact, 
the Sifting Process can recover low amplitude riding waves, 
which may appear as gentle humps in the original signal, 
65 with repeated siftings. 

Still another complication is that the envelope mean 40 
may be different from the true local mean for nonlinear data. 
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Consequently, some asymmetric wave forms can still exist 
no matter how many times the data are sifted. This must be 
accepted because, after all, the inventive method is an 
approximation as discussed before. 

Other than these theoretical difficulties, on the practical 
side, serious problems of the spline fitting can occur near the 
ends, where the cubic spline being fit can have large swings. 
Left by themselves, the end swings can eventually propagate 
inward and corrupt the whole data span especially in the low 
frequency components. A numerical method has been 
devised to eliminate the end effects details of which will be 
given later. Even with these problems, the Sifting Process 
can still extract the essential scales from the data. 

The Sifting Process serves two purposes: to eliminate 
riding waves and to make the wave profiles more symmetric. 
Toward these ends, the Sifting Process has to be repeated. 
Because only the first component signal tq has been gener- 
ated so far, the decision step 260, which tests successive 
component signals to see if they satisfy the definition of an 
IMF, is bypassed during the first iteration. 

Thus, step 265 is performed which treats the component 
signal as the physical signal in the next iteration. The next 
iteration is then performed by executing steps 200-250. In 
step 250, the second component signal h ±1 is generated by 
subtracting the envelope mean from the physical signal (in 
this iteration, the first component signal h 2 is treated as the 
physical signal). In more formal terms: 
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the resulting IMF a pure frequency modulated signal of 
constant amplitude. 

To guarantee that the IMF component retains enough 
physical sense of both amplitude and frequency 
5 modulations, a stopping criterion is employed to stop the 
generation of the next IMF component. This stopping cri- 
terion is part of the computer implemented method and is 
shown as step 260 in FIG. 3 (b). Step 260 is a decision step 
that decides whether the stopping criteria has been satisfied. 
The preferred stopping criteria determines whether three 
successive component signals satisfy the definition of IMF. 
If three successive component signals all satisfy the defini- 
tion of IMF, then the Sifting Process has arrived at an IMF 
and should be stopped by proceeding to step 270. If not, step 
15 260 starts another iteration by proceeding to step 265 as 
described above. 

Alternatively, another stopping criteria could be used that 
determines whether successive component signals are sub- 
stantially equal. If successive component signals are sub- 
2Q stantially equal, then the Sifting Process has arrived at an 
IMF and should be stopped by proceeding to step 270. If not, 
step 260 starts another iteration by proceeding to step 265 as 
described above. 

Determining whether successive component signals are 
25 substantially equal in the alternative stopping criteria limits 
the size of the standard deviation, sd, computed from the two 
consecutive sifting results as follows: 


FIG. 4(h) shows the second component signal h 1± . 
Although the second sifting shows great improvement in the 
signal with respect to the first sifting, there is still a local 
maximum below the zero line. After a third sifting, the result 
(third component signal h 12 ) is shown in FIG. 4(h). Now all 
the local maxima are positive, and all the local minima are 
negative, but to ensure this configuration is stable, the 
Sifting Process should be further repeated. In general, the 
Sifting Process is repeated at least 3 more times and, in 
general, K times to produce h lk . If no more new extrema are 
generated, then h lk , is an IMF. In formal terms: 

The resulting first IMF component is shown in FIG. 4(f) 
after 9 siftings. The first IMF component of the physical 
signal may be designated as such in step 270 and stored in 
step 275 in memory 70, c-^h-^. 

As mentioned above, all these manipulations are carried 
out numerically in a computer 52. There is not explicit close 
form analytic expression for any of the computer imple- 
mented steps. 

As described above, the process is indeed like sifting of 
the data by the computer 52 because it separates the finest 
(shortest time scale) local mode from the data first based 
only on the characteristic time scale. The Sifting Process, 
however, has two effects: 

a. to eliminate riding waves, and 

b. to ensure the envelopes generated by maxima and 
minima are symmetric. 

While the first condition is necessary for the instantaneous 
frequency to be meaningful (as discussed below), the second 
condition is also necessary in case the neighboring wave 
amplitudes have too large a disparity. Unfortunately, the 
effects of the second condition, when carried to the extreme, 
could obliterate the physically meaningful amplitude fluc- 
tuations. Therefore, the Sifting Process should be applied 
with care, for carrying the process to an extreme could make 


30 


sd 


■2 


I ^l(it-l)(0 - hlk(t)) I 

A ?(*- 1)(0 


Avery rigorous and preferred value for sd is set between 0.2 
and 0.3. Of course, if faster processing is desired, then a 
35 trade -off such as a less rigorous value for sd may be used. 

Overall, the first IMF component c 1 should contain the 
finest scale or the shortest period component of the physical 
signal 10. 

Before extracting the next IMF component, a test should 
40 be conducted to determine if the Sifting Process should stop. 
The stopping criteria is shown in Step 300. Step 300 
determines whether the component signal has more than 2 
extrema. If not, all of the IMF’s have been extracted and the 
Sifting Process is stopped by proceeding to step 310. If so, 
45 then additional IMF’s may be extracted by continuing the 
process in step 320. 

Step 270 recognizes that an IMF component has been 
successfully generated by the Sifting Process by setting the 
component signal equal to an intrinsic mode function. More 
50 specifically, step 270 causes the computer 52 to store the 
component signal h lk as an intrinsic mode function in 
memory 70. 

Then, the first IMF is Separated from the physical signal 
in step 290 to generate a residual signal. In particular, a 
55 residual signal is generated by subtracting the intrinsic mode 
function from the physical signal. In formal terms: 

m-ci= ri . 

Because the residue, r 1? still includes information of longer 
60 period components, it is treated as the new physical data and 
subjected to the same Sifting Process as described above. 
Step 320 performs this function by treating the residual 
signal as the physical signal in the next iteration. Thereafter, 
the next iteration is performed beginning with the execution 
65 of step 200 as described above. 

The Sifting Process is repeated for all the Subsequent r^’s. 
This iterative procedure may be expressed as: 
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ri- c 2 =r 2 , 

r n-l~ C n =r n‘ 

In step 300 the Sifting Process ends by proceeding to stop 
step 310 if the residual signal r M has more than 2 extrema. 5 
Otherwise, the method proceeds to step 320. 

In other words, Step 310 stops the Sifting Process if the 
residual signal r n is monotonically increasing or decreasing. 
This stopping criterion is based on the fact that an IMF 
cannot be extracted from a mono tonic function. Even for 10 
data with zero mean, the final residue still can be different 
from zero. For data with a trend, the final residue should be 
that trend. If r n is not monotonically increasing/decreasing, 
then a next iteration is performed beginning with step 320. 

In summary, the Sifting Process decomposes the physical 15 
signal X(t) into a series of intrinsic mode functions and a 
residue which may be expressed as: 

x(t) = c,- + r n . 20 


In other words, the invention extracts a series of IMFs by 
Sifting the physical signal with a computer implemented ^ 
Empirical Mode Decomposition method. This IMF series 
cannot be generated or derived by any analytic method. It 
can only be extracted by the invention through a specially 
programmed computer through the Sifting Process (local 
extrema or curvature extrema type). 

Display of Selected Results 

The various results of the above-described computer 
implemented method are displayed in step 116 . These dis- 
plays are extremely useful in analyzing the underlying 
physics of the physical phenomenon being studied as 
described above. Furthermore, particular examples of these 
displays and the increased understanding of the underlying 
physics which these displays permit are discussed in the 
following section. 

For example, the inventions generates various Hilbert 
spectra displays in the display step 116 . As mentioned 
above, both color coded maps and contour maps may be 
employed to display the Hilbert spectra in display step 116 . 

In addition, the color coded maps convey information to the 
operator in a uniquely accessible way permitting a more ^ 
thorough and deeper understanding of the physical phenom- 
enon and may be considered as necessary to analyze some 
physical phenomena. 

The displays generated by the invention in display step 
116 are central to the invention because they allow an 
operator to analyze the underlying physics of the physical 
phenomenon being studied. 

The display step 116 outputs displays to display 62 . As 
mentioned above, display 62 includes devices such as a 
cathode ray tube and a flat panel display. As an alternative, 
display step 116 may generate a hard copy output by 
utilizing printer 64 or may send the generated display to 
output device 66. 

Application of the Preferred Embodiment 

60 

Here, for ease in discerning an Antarctic circumpolar 
wave in the perimeter of the ice pack, a time series of the sea 
ice extents (the sum of the areas of pixels with ice concen- 
trations of 15% or more) is constructed in 1 -degree longi- 
tudinal sectors for the period 1978-1996, as observed with 65 
multichannel microwave imagers on board the NASA Nim- 
bus 7 and the DOD DMSP F8, Fll and F13 satellites the 
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data are transmitted to input device 58. An example of the 
data is shown in FIG. 5 for a typical day in Austral Winter. 

Next, Complex Singular- Value Decomposition (CSVD) is 
used to decompose the time series of the 360 sectors into its 
principal components, after subtracting its linear trend line, 
effectively separating the spatial and temporal values. The 
results of the CSVD procedure can be expressed as: 

y(3322, 360)= £7(3322,3 60) * 5(360,360) * V(360,360)' 

after the original data have been subjected to a Hilbert 
transform to produce a complex data vector matrix, Y 
referred to as A complex hereinabove. Y is represented by the 
product of two complex unitary matrices, U, which contains 
the temporal information for each of the 360 CPCs (PC#) 
and V, which contains the spatial information for each of the 
CPCs. S is a diagonal matrix of real eigenvalues. The prime 
on V denotes the matrix transpose. The time number 3322 
corresponds to the number of data days in the 18.2-year time 
series (including a 21 -day data day gap in the DMSP F8 
imager record, which has been filled by interpolation with a 
modeled seasonal cycle). 

The spatial and temporal factors of the real part of the first 
three principal components are shown in FIGS. 6(a)-(b). 
Plots of the eigenvalue ratios of the CPCs and also of the 
cumulative fraction of the signal in the first n modes 
(obtained by taking the ratio of the sum of the eigenvalues 
of the first n modes and the sum of all the eigenvalues) are 
shown in FIGS. 6(a) and 6(b). FIG. 6(b) indicates that about 
34% of the signal is contained in the first 3 CPCs. FIG. 6(a) 
indicates that the amplitude of the fourth CPC is about 11% 
of the first. Thus, we need utilize only the first three CPCs 
to obtain a reasonable approximation to the original signal 
and to visualize any circumpolar wave that might be present. 

The temporal portions of the first 3 CPCs are processed by 
means of Empirical Mode Decomposition (EMD) into their 
intrinsic modes (IM), resulting in a collection of 2-3 CPCs 
for each intrinsic mode. It is important to note that the 
presence of several intrinsic modes for each CPC indicates 
that they are multi-modal. In particular, visual inspection of 
the temporal factor of CPC1 (FIG. 9(b)) might lead one to 
conclude that it contains only the seasonal cycle, but as is 
shown later this is far from the case. A simplistic description 
of this method of sifting a data vector, y(n), into nearly 
orthogonal intrinsic modes follows: 

1. Obtain two vectors containing the extrema of y(n). 

2. Fit maxima and minima with two cubic spline vectors of 
length n. 

3. m 1 (n)=average of the two cubic spline vectors. 

4. First component, c 1 =y(n)-m 1 (n). 

5. Repeat steps 1-3 on m 1 (n) to obtain m 2 (n). 

6. Second component c 2 =m 1 (n)-m 2 (n). 

i. ith component, c~m,_ 1 (n)-m 2 (n). 

Continue until c, contains less than one whole oscillation. 
This usually converges for i=8 to 10. 

Conditional tests are included for assuring that the num- 
ber of zero -crossings equal the number of extrema, for how 
many time series points are allowed between the extrema of 
the various modes, or how many are allowed between the 
extrema of the curvatures of the modes. By selecting input 
criteria judiciously the resultant modes are nearly orthogo- 
nal. A benefit of this approach is that the oscillations in the 
various modes are confined to narrow instantaneous fre- 
quency bands, greatly simplifying the interpretation of the 
results. Another benefit is that the original signal (within 
machine accuracy) is regained when the modes are added 
back together, indicating that spurious signals have not been 
introduced inadvertently as a result of the decomposition. 
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The intrinsic modes of CPC1 are shown in FIG. 7. CPC1 
appears to have 8 distinct oscillatory modes. The first mode 
contains the day-to-day fluctuations in the sea ice extent 
attributable to fluctuations in the weather. There is a flat area 
surrounded by two transients in this record occurring at the 
time of the aforementioned 21 day data gap, reflecting the 
lack of single-day fluctuations in the interpolated data. The 
second mode is largely semi-monthly oscillations, possibly 
related to the 14-day tidal frequency. The third mode appears 
to be largely monthly oscillations, with an irregular ampli- 
tude modulation envelope of an annual cycle. These oscil- 
lations were interpreted as the influence of the monthly 
oceanic tides on the edge of the sea ice, and the minima of 
the envelope reflecting the reduction of signal when the ice 
is at its annual minimum extent. The fourth mode is a 
mixture of periods larger than one month and shorter than 
one year of unknown origin. The fifth mode is what we shall 
call the seasonal, rather than annual, cycle because it is not 
a monochromatic sinusoid. The shape of the fifth mode 
differs from other seasonal cycle models, where 5 harmonics 
of the annual cycle are weighted and summed. The shape of 
the fifth mode derived using the present invention was found 
to approximate closely the observed seasonal cycles. The 
difference is attributable to the inability of Fourier analysis 
to produce a meaningful spectrum for a non-stationary 
phenomenon. 

Modes 6-8 are of the most interest in the context of a 
search for a circumpolar wave in the Antarctic sea ice 
perimeter. The amplitudes of modes 6-8 are (respectively) 
about 0.25, 0.1, and 0.05 of the seasonal cycle amplitude 
(Mode 5). Modes 6 and 7 are quasi-biennial and quasi- 
quadrennial periods similar to those observed earlier in the 
Antarctic sea ice canopy with band-limited Fourier analysis 
of the 9-year Nimbus 7 data record and attributed to an 
atmospheric teleconnection with the El Nino and Southern 
Oscillation phenomena. Mode 8 contains periods in the 
range 6-7 years, which was only hinted at in the earlier 
studies of a 9 -year data set. Finally, the last mode (9) 
contains the residual trend of the data set remaining after the 
initial subtraction of the trend line obtained by ordinary least 
squares. Although it displays some curvature. Mode 9 can- 
not be described convincingly as oscillatory. 

The similarity between the modal structure of each of the 
first 3 CPCs and the rapid fall-off of their eigenvalues permit 
reconstruction of the original signal from its first 3 CPCs to 
a good approximation. Furthermore, the intrinsic modes of 
the original data can be reconstructed to a good approxima- 
tion from the individual intrinsic modes of the first 3 CPCs. 

As noted above, the data were reconstructed in two 
different ways. First, to simulate the 3-7 year bandpass filter 
and added together the intrinsic modes 7-9 of CPC1 (and 
similar modes of CPC2 and 3) to reconstruct a low-pass 
filtered time series. This excludes the quasibiennial modes 
(Mode 6 of CPC1) as well as the seasonal and shorter-period 
modes. Second, the quasiquadrennial (QQ) mode (Mode 7 
of CPC1) is isolated to look for differences between the 
QQ -filtered and lp -filtered results. 

Time/longitude contours of the lp-filtered data (FIG. 9(a)) 
show eastward-propagating azimuthal motion, but no com- 
plete circumpolar circuit. The most prominent one begins in 
about 1990, with the crest at about 220° E and trough at 
about 320° E. The crest propagates eastward until about 
1994, when it reaches about 40° E. The trough reaches 40° 
E in about 1992. The complicated combination of stationary 
and propagating troughs may be seen from a different 
perspective in the 3-dimensional contour diagram in FIG. 
9(b). There seems to be a barrier between about 120-150° E 
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(the part of the Southern Ocean south of the western Pacific) 
that the propagating waves do not cross. 

The aforementioned barrier appears to be absent in the 
propagation patterns of the QQ modes shown in FIGS. 
5 10(a)-(Z?). Focusing again on two of the most prominent 
features, a trough beginning in about 1984 at about 220° E 
propagates eastward to 360° E until about mid-1986, then 
regresses westward for about Vi a year, resuming the east- 
ward propagation to reach 360° E once again in about 1994, 
10 finally stopping at about 20° E in about mid-1995. Thus, the 
trough propagated for about 1.5 circumpolar cycles at a rate 
of about 1 cycle in 8 years. The adjacent crest also starts at 
about 220° E but in about mid- 1985. It propagates eastward 
crossing the Greenwich Meridian in about 1985. There is a 
15 brief westward regression in about mid-1989, after which it 
continues eastward to about 300° E in about mid-1995 at 
which point the crest begins another westward regression. 

Eastward propagation of a circumpolar wave have been 
reported in sea ice extents as well as in sea level pressures, 
20 and sea surface winds and temperatures. The data were 
preprocessed by passing through a band-pass filter with 
admittance from 3 to 7 years. Sea ice extents were analyzed 
in 5° meridional sectors, and other parameters along the 56° 
S parallel. The contour patterns were interpreted on their 
25 resulting position, time diagrams as circumpolar waves 
propagating eastwards and completing the circular traverse 
in about a decade. Thus, the low-pass filter procedure is 
believed to encompass the entire band-passed filtered oscil- 
latory behavior. 

30 Initially, this Antarctic sea ice extent investigation was 
undertaken to confirm what was the commonly held belief 
and to examine in detail what has been previously reported, 
regarding the circumpolar wave in sea ice extents. 
Unexpectedly, however, the sea ice extent results obtained 
35 with lp filter of the present invention do not agree with prior 
results. Perhaps, the prior practice of binning of the sea ice 
extents into 5° sectors, thereby limiting the spatial resolution 
to about 10° of longitude, may have bridged the “barrier” 
between 120-150° that was observed in lp-filtered results of 
40 the present invention. Alternatively, the high frequency end 
of the prior art band-pass filter with a half-power point of 3 
years, may suppress part of the signal admitted by lp-filter 
of the present invention. When oscillations of 3 years or less 
are eliminated with QQ filter, the results are much more in 
45 agreement. 

The results of applying the present invention to this 3D 
physical phenomenon clearly demonstrate the necessity of 
using a complex number representation of the time series in 
order to retain the circumpolar motion in the sea-ice pack 
50 when applying singular value decomposition. Earlier 
attempts at observing the circumpolar motion with real 
numbers was unsuccessful. Accordingly, it is clear that care 
must be used in selecting the appropriate filters and the need 
for high spatial resolution in order not to lose important 
55 details in the wave propagation. Finally, it is imprudent to 
presume that complex principal components in themselves 
represent fundamental modes of an oscillatory system. For 
instance, the first CPC shown here in FIGS. 6(a)-(b) could 
easily be mistaken by inspection for the season cycle. While 
60 the seasonal cycle is the largest oscillatory component of the 
first CPC, the first CPC also has 7 other significant oscilla- 
tory modes. 

Alternative Embodiments 

As described above, in the EMD step 110, upper and 
65 lower envelopes 20, 30 are constructed with a cubic spline 
in steps 210 and 230, respectively. This cubic spline fitting, 
however, has both overshoot and undershoot problems. 
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These problems can be alleviated by using sore sophisticated 
spline methods, such as the taut spline in which the tension 
of the spline curve can be adjusted. 

Another alternative is higher-order spline fitting. 
Although such higher-order spline fitting may be more 
accurate, it will, however, introduce more inflection points 
or extrema, and consume more computing time. Therefore, 
it is not recommended as a standard operation. Only in 
special cases, it may be tried. 

As the spline fitting procedure is time consuming, more 
efficient methods can be devised by using simple mean 
values of successive extrema instead of computing the 
envelope-mean. In this way, only one spline fitting is 
required rather than two. Although this alternative is easier 
and faster to implement, the shortcomings are more severe 
amplitude averaging effects when the neighboring extrema 
are of different magnitudes. The successive-mean method 
will have a stronger forcing to reach uniform amplitudes, in 
which the true physics associated with amplitude will be 
destroyed. Therefore, the successive-mean method should 
only be applied where the amplitudes of the physical signal 
components are constants. 

Either the envelope mean or the successive mean method, 
when applied with the requirement of absolute symmetry, 
will produce the absurd result of uniform amplitude IMF’s. 
Therefore, the criteria in the Sifting Process should be 
chosen judiciously. One should avoid too stringent a crite- 
rion that we would get uniform amplitude IMF’s. On the 
other hand, one should also avoid too loose a criterion that 
would produce components deviating too much from IMF’s. 

It is well known that the most serious problem of spline 
fitting is at the ends, where cubic splines can have wide 
swings if left unattended. As an alternative, the invention 
may utilize a method of adding characteristic waves at the 
ends of the data span. This confines the large swings 
successfully. 

The method of adding characteristic waves is not con- 
ventional. In contrast, the conventional window that is often 
applied to Fourier transform data results in loss of useful 
data. To avoid this data loss and to confine swings at the ends 
of the data span, the invention extends the data beyond the 
actual data span by adding three additional characteristic 
waves at each end of the data span. 

The characteristic waves are defined by the last wave 
within the data span at the end of the data span. In other 
words, a characteristic wave is added to each end of the data 
span having an amplitude and period matching the last wave 
within the data span. This characteristic wave is a sinusoidal 
waveform that is extended three sinusoidal wave periods 
beyond the data span at each end. This process is repeated 
at the other end of the data span. In this way, spline fitting 
at the end of the data span, which can otherwise have a wide 
swing, is confined. In other words, by adding the extra 
characteristic waves at the ends beyond the data span, the 
spline curve will be tied down so that it will not have wild 
or excessive swings that would otherwise corrupt the data 
processing and analysis that utilizes these cubic splines. 

Other than the spline fitting, the Hilbert transform may 
also have end effects. Because the first and the last points of 
the data are usually of different values, the Fourier transform 
will introduce additional components to bridge over the 
difference resulting in the well-known Gibbs phenomena. To 
eliminate it in the Fourier transform, various windows have 
been adopted (see, for example, Brigham, 1974, “The fast 
Fourier Transform ”, Prentice-Hall, Englewood Cliffs, N.J.). 

Instead of a window which will eliminate some useful 
data at the end, the invention again adds two characteristic 
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waves at either end. These waves all start from zero at the 
beginning, and end at zero at the end. Thus, the annoying 
Gibbs phenomena are greatly reduced. 

Still further, the Hilbert transform needs over-sampled 
5 data to define the instantaneous frequency precisely. In 
Fourier analysis, the Nyquist frequency is defined by two 
points per wave, but the frequency is defined for a wave 
covering the whole data span. In the invention, the instan- 
taneous frequency is defined through a differentiation 
10 process, and thus more data points will help defining the 
frequency more accurately. Based on the inventor’s 
experience, a minimum number of data points to define a 
frequency is five (or four At’s). The lack of fine time steps 
can be alleviated through interpolating more points between 
15 available data. As a spline interpretation would also not 
create nor annihilate scales, it can also be used for the 
interpolation when the data are very jagged from under- 
sampled data. The smoothed data though have a longer 
length and are sometimes easier to process. The interpola- 
20 tion may give better frequency definition. 

Particular Limitations of The Invention 

The dependence on the existence of scale for mode 
definition has one limitation: the decomposition method 
cannot separate signals when their frequencies are too close. 
25 In this case, there would not be any characteristic scale: 
therefore, physically they are identical. This may be the most 
severe limitation of the invention, but even here, the 
invented method can still work as well as the Fourier 
Analysis. 

30 Particular Advantages of The Invention 

The strength of the EMD method should be reiterated. 
EMD is built on the idea of identifying the various scales in 
the data which are quantities of great physical significance. 
Therefore, in the local extrema and curvature extrema Sift- 
35 ing Processes, orthogonality is not a consideration, but 
scales are. Since orthogonal decomposition is a character- 
istic for linear systems, violating this restriction is not a 
shortcoming but a breakthrough. Therefore, the decomposed 
IMF’s may or may not be orthogonal. As such, this method 
40 can be applied to nonlinear data. Though the IMF’s in most 
cases are practically orthogonal, it is a coincidence rather 
than a requirement of the EMD. 

Finally, though the EMD method will give IMF 
components, the individual component does not guarantee 
45 well-defined physical meaning. This is true for all 
decompositions, especially for the methods with a priori 
basis. In most cases, however, the IMF’s do carry physical 
significance. Great caution should be exercised in making 
such attempts. The rule for interpreting the physical signifi- 
50 cance of the IMF’s is that the scales should be clearly 
separated. Together with the Hilbert spectrum, the totality of 
the presentation should give a much more detailed repre- 
sentation of the physical processes than conventional meth- 
ods. 

55 The invention being thus described, it will be obvious that 
the same may be varied in many ways. Such variations are 
not to be regarded as a departure from the spirit and scope 
of the invention, and all such modifications as would be 
obvious to one skilled in the art are intended to be included 
60 within the scope of the following claims. 

I claim: 

1. A computer implemented method of analyzing three- 
dimensional (3D) phenomena, said method comprising the 
steps of: 

65 inputting a representation of a 3D phenomenon; 

creating a complex representation of said 3D phenom- 
enon; 
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recursively sifting said complex representation using 
Empirical Mode Decomposition to extract an intrinsic 
mode function indicative of an intrinsic oscillatory 
mode; and 

processing the intrinsic mode function. 5 

2. A computer implemented method of analyzing 3D 
phenomena as in claim 1, wherein the step of creating a 
complex function comprises the steps of: 

reshaping said 3D representation into a two-dimensional 
(2D) representation; and 10 

applying a Hilbert transform to said reshaped 3D 
representation, the result being said complex represen- 
tation. 

3. A computer implemented method of analyzing 3D 
phenomena as in claim 2, wherein the analyzed 3D phe- 
nomenon is a seasonal physical phenomenon and the step of 
reshaping the 3D representation further comprises: 

deseasoning said reshaped representation (AA), the Hil- 
bert transform being applied to said deseasoned 
reshaped representation (AA Deseasoned ). 

4. A computer implemented method of analyzing 3D 
phenomena as in claim 3, wherein the step of deseasoning 
the 2D representation further comprises removing a trend 
line using multiple linear regression, the deseasoned 2D 
representation being described by the relationship: 

25 

5 

AA Deseasoned = AA - a 0 + a x t + ^ [a 2k cos{27rkt / t) + 0^+1 )Sm(2nkt / t)], 
k= 1 

wherein a 0 , a 1? a 2Ar , a 2/r+1 result from multiple linear regres- 30 
sion and x is the annual period. 

5. A computer implemented method of analyzing 3D 

phenomena as in claim 3, wherein in the reshaping step the 
3D representation is restructured as a series of longitudinal 
sector extents. 35 

6. A computer implemented method of analyzing 3D 

phenomena as in claim 5, wherein said series is a time series 
and wherein each of said longitudinal extents is the sum of 
all areas of a selected granularity having a selected charac- 
teristic identified as being above a given threshold. 40 

7. A computer implemented method of analyzing 3D 
phenomena as in claim 6, wherein said selected area of 
granularity is represented by a pixel on a display of said 
phenomenon and said threshold is 15%. 

8. A computer implemented method of analyzing 3D 45 
phenomena as in claim 3, wherein the step of applying the 
Hilbert transform preserves spatial phase information in said 
2D representation. 

9. A computer implemented method of analyzing 3D 
phenomena as in claim 8, said method further comprising: 

applying singular value decomposition to the result of the 50 
Hilbert transform, complex principal components 
being extracted from the Hilbert transform results. 

10. A computer implemented method of analyzing 3D 
phenomena as in claim 9, wherein the result of the Hilbert 
transform may be represented as a two-dimensional matrix, 55 
A complex* the extracted spatial component may be repre- 
sented as a two-dimensional matrix, U, and the result of the 
extracted temporal component can be represented as a 
two-dimensional matrix, V, and wherein 

60 

A com p/&c=U * S * V' , 

S being an eigenvalue matrix and V' being the transpose of 
V. 

11. A computer implemented method of analyzing 3D 
phenomena as in claim 10, wherein EMD is performed on 65 
the temporal parts (V) of a plurality of complex principal 
components. 
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12. A computer implemented method of analyzing 3D 
phenomena as in claim 11, wherein the plurality of complex 
principal components is 10. 

13. A computer implemented method of analyzing 3D 
phenomena as in claim 12, wherein, in the processing step, 
the 3D phenomena representation is reconstructed by select- 
ing one or more modes, selected said one or more modes 
constituting a band-pass filter. 

14. A computer implemented method of analyzing 3D 
phenomena as in claim 13, wherein the reconstructed filtered 
(A complex , filtered) 3D phenomenon may be represented as: 

Complex, fii t erecrU(position, 1 : 10 ) * 5 ( 1 : 10 , l:10)*V filtere Jtime, 

1 : 10 )'. 

15. A computer implemented method of analyzing 3D 
phenomena as in claim 13, wherein the real part of the 
reconstructed representation is reshaped into a sequence of 
filtered grid maps. 

16. A computer implemented method of analyzing 3D 
phenomena as in claim 15, wherein the processing step 
further comprises displaying filtered grid maps. 

17. A three-dimensional (3D) phenomena analysis appa- 
ratus comprising: 

a receiver receiving a representation of a 3D phenom- 
enon; 

a component extractor extracting a temporal component 
and a spatial component to form a complex represen- 
tation of said 3D phenomenon; 

an intrinsic mode extractor recursively sifting said com- 
plex representation using Empirial Mode Decomposi- 
tion to extract an intrinsic mode function indicative of 
an intrinsic oscillatory mode; and 

means for processing the intrinsic mode function to 
regenerate a displayable representation of said 3D 
phenomena. 

18. A 3D phenomena analysis apparatus as in claim 17, 
wherein the component extractor comprises: 

a reshaper reshaping said 3D representation into a two- 
dimensional (2D) representation; and 

a Hilbert Spectrum generator applying a Hilbert transform 
to said 2D representation, the result being said complex 
representation. 

19. A 3D phenomena analysis apparatus as in claim 18, 
wherein the reshaper comprises: 

a segmenter segmenting the 3D representation into a 
plurality of uniform segments, the area of each of said 
segments being a plurality of uniform granularity areas; 
and 

a counter counting selected ones of said uniform granu- 
larity areas. 

20. A 3D phenomena analysis apparatus as in claim 19, 
wherein the reshaper further comprises: 

weight determination means for determining a weight 
attached to each of said uniform granularity areas, said 
counter counting said uniform granularity areas having 
a weight above a selected threshold. 

21. A 3D phenomena analysis apparatus as in claim 18, 
wherein the analyzed 3D phenomenon is a seasonal physical 
phenomenon and the component extractor further com- 
prises: 

a deaseasoner deseasoning said reshaped representation, 
said deseasoned reshaped representation being said 2D 
representation (AA). 

22. A 3D phenomena analysis apparatus as in claim 21, 
wherein the deseasoner comprises: 
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sinusoidal calculation means for determining time depen- computer readable program code means for receiving a 

dent sine and cosine values; representation of a 3D phenomenon; 

a multiplier combining selected coefficient with calcu- computer readable program code means for forming a 

lated sine and cosine values: and complex representation of said 3D phenomenon, said 


an adder combining a plurality of products from said 
multiplier with said reshaped representation. 

23. A 3D phenomena analysis apparatus as in claim 21, 
wherein said deseasoner includes multiple linear regression 
means and also removes trend lines from said reshaped 
representation (AA), said deseasoned reshaped representa- 
tion (AA Deseasoned ) having the form: 

5 

A A Deseasoned = AA - dp + djt + [a 2k cos(27rkt / t) + d {2k+l) sm(27rkt / r)], 

Jt=l 

multiple linear regression means generating a 0 , a 1? a 2Ar and 
a 2 *+i an d? T being the annual period of said 3D representa- 
tion. 

24. A 3D phenomena analysis apparatus as in claim 21, 
wherein the Hilbert spectrum generator preserves spatial 
phase information in said 2D representation. 

25. A 3D phenomena analysis apparatus as in claim 24, 
further comprising: 

a singular value decomposer extracting complex principal 
components from the complex representation. 

26. A 3D phenomena analysis apparatus as in claim 25, 
wherein the Hilbert spectrum generator generates a two- 
dimensional matrix, A compIex , the singular value decomposer 
extracts a spatial component represented as a two- 
dimensional matrix, U, and a temporal component repre- 
sented as a two-dimensional matrix, V, and wherein 

A cowp/e *=u*s*v, 

S being an eigenvalue matrix and V' being the transpose of 

V. 

27. A 3D phenomena analysis apparatus as in claim 26, 
wherein the intrinsic mode extractor performs EMD on the 
temporal parts (V) of a plurality of complex principal 
components. 

28. A 3D phenomena analysis apparatus as in claim 27, 
wherein the processing means selects one or more modes to 
reconstruct the filtered 3D phenomena representation. 

29. A 3D phenomena analysis apparatus as in claim 28, 

wherein the reconstructed filtered (A complex 3D phe- 

nomenon is represented as 

Acomplex, fiitered=U(position, l:f)*V filtered (time, 1 :f)\ 

f being the number of complex principal components 
included in the filter. 

30. A 3D phenomena analysis apparatus as in claim 28, 
wherein the processing means reshapes the real part of the 
reconstructed representation into a sequence of filtered grid 
maps. 

31. A 3D phenomena analysis apparatus as in claim 30, 
further comprising: 

a display displaying filtered grid maps. 

32. A 3D phenomena analysis apparatus as in claim 17, 
further comprising: 

a display displaying said regenerated representation of 
said 3D phenomenon. 

33. A computer program product for analyzing 3D 
phenomena, said computer program product comprising a 
computer usable medium having computer readable pro- 
gram code comprising: 


5 complex representation including a temporal compo- 
nent and a spatial component; 

computer readable program code means for recursively 
sifting said complex representation using Empirical 
Mode Decomposition to extract an intrinsic mode func- 
10 tion indicative of an intrinsic oscillatory mode; and 

computer readable program code means for processing 
the intrinsic mode function to regenerate a displayable 
representation of said 3D phenomena. 

34. A computer program product for analyzing 3D phe- 
5 nomena as in claim 33, wherein the computer readable 

program code means for forming the complex representation 
comprises: 

computer readable program code means for reshaping 
said 3D representation into a two-dimensional (2D) 
representation; and 
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computer readable program code means for applying a 
Hilbert transform to said 2D representation, the result 
being said complex representation. 

35. A computer program product as in claim 34, wherein 
the computer readable program code means for reshaping 
comprises: 

computer readable code means for restructuring said 3D 
representation as a series of longitudinal extents. 

36. A computer program product as in claim 35, wherein 
30 each of the longitudinal extents represents a proportion of a 

sector of said 3D representation, said proportion indicating 
a portion of said sector having a measured property beyond 
a threshold level. 

37. A computer program product for analyzing a 3D 
35 phenomena as in claim 34, wherein the analyzed 3D phe- 
nomenon is a seasonal physical phenomenon and the com- 
puter readable program code means for extracting compo- 
nents further comprises: 

computer readable program code means for deseasoning 
40 said reshaped representation, said deseasoned reshaped 
representation being said 2D representation (AA). 

38. A computer program product for analyzing a 3D 
phenomena as in claim 37, wherein the computer readable 
program code means for deseasoning the 2D representation 

45 removes trend lines and comprises computer readable pro- 
gram code means for applying multiple linear regression, the 
deseasoned 2D representation (AA Deseasone J) satisfying the 
relationship: 

50 i 

AA Deseasoned = AA-a 0 +a l t+ 2^ [d 2k cos{2nkt j t) + a {2k+l) sm{2jrkt / t)\, 

k = l 

a 0 , a 1? a 2 £ and a 2 ^ +1 being generated by multiple linear 
55 regression and, x being the annual period of said 3D repre- 
sentation. 

39. A computer program product for analyzing 3D phe- 
nomena as in claim 38, wherein the computer readable 
program code means for applying the Hilbert transform 

60 preserves spatial phase information in said 2D representa- 
tion. 

40. A computer program product for analyzing 3D phe- 
nomena as in claim 39, further comprising: 

computer readable program code means for applying 
65 singular value decomposition to the result of the Hilbert 
transform, complex principal components being 
extracted from the Hilbert transform results. 
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41. A computer program product for analyzing 3D phe- 
nomena as in claim 40, wherein the result of the Hilbert 
transform may be represented as a two-dimensional matrix, 

A complex > the extracted spatial component may be repre- 
sented as a two-dimensional matrix, U, and the result of the 5 
extracted temporal component can be represented as a 
two-dimensional matrix, V, and wherein 

^complex= U*S*V', 

S being an eigenvalue matrix and V' being the transpose of 
V. 

42. A computer program product for analyzing 3D phe- 
nomena as in claim 41, wherein the computer program 
product means for recursively sifting applies EMD to the 
temporal parts (V) of a plurality of complex principal 
components. 

43. A computer program product for analyzing 3D phe- 
nomena as in claim 42, where the computer readable pro- 
gram code means for the processing further comprises: 2Q 
computer program product means for reconstructing the 3D 
phenomena representation by selecting one or more modes, 
selected said one or more modes constituting a band-pass 
filter. 
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44. A computer program product for analyzing 3D phe- 
nomena as in claim 43, where the computer readable pro- 
gram code means for reconstructing reconstructs the filtered 
3D phenomenon (A compIeXi fiIlered ) represented as: 

A complex, fiitered = U (position, l:f)*S(l:f, 1 :f ) * V fi lter erf (time, 1 

f being the number of complex principal components 
included in the filter. 

45. A computer program product for analyzing 3D phe- 
nomena as in claim 43, wherein the computer readable 
program code means for processing comprises computer 
program product means for reshaping the real part of the 
reconstructed representation into a sequence of filtered grid 
maps. 

46. A computer program product for analyzing 3D phe- 
nomena as in claim 45, wherein the computer readable 
program code means for processing further comprises com- 
puter readable program code means for providing filtered 
grid maps for display. 



